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Abstract 

A complete method is proposed to compute a certified, or ambient isotopic, meshing for 
an implicit algebraic surface with singularities. By certified, we mean a meshing with correct 
topology and any given geometric precision. We propose a symbolic-numeric method to com- 
pute a certified meshing for the surface inside a box containing singularities and use a modified 
Plantinga-Vegter marching cube method to compute a certified meshing for the surface inside a 
box without singularities. Nontrivial examples are given to show the effectiveness of the algo- 
rithm (see Fig. [T]!. To our knowledge, this is the first method to compute a certified meshing for 
surfaces with singularities. 

Keywords. Surface, curve, topology, ambient isotopic meshing, marching cube, symbolic com- 
putation, interval arithmetic. 



1 Introduction 




Figure 1 : Isotopic meshing for surfaces with singular points and singular curves 

To determine the topology of a given algebraic surface and to use triangular meshes to approxi- 
mately represent the surface are fundamental operations in computer graphics and geometric model- 
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ing. Meshing of surfaces could be used to display the surface correctly and to perform engineering 
applications on the surface, such as the finite element analysis. A survey on this topic can be found 
in a. 

We consider an implicit surface defined by /(x, y,z) =0 where f{x, y, z) is a square free polyno- 
mial with rational numbers as coefficients. There exists a large amount of work on meshing implicit 
surfaces. Please see the work |fTl |4l |26l and the literatures cited in them. Recent work focuses on 
isotopic meshing [5|. Simply speaking, a meshing is called isotopic if it has the same topology and 
the same geometry as the surface (for definition see Section 2). A meshing is called ambient isotopic 
or certified if it is isotopic and approximates the surface to any give precision. There exist four main 
approaches to compute isotopic meshings for surfaces: the marching cube method, the Morse the- 
ory method, the Delaunay refinement method, and the CAD (Cylindrical Algebraic Decomposition) 
based method. 

The famous marching cube method repeatedly subdivides the space into smaller cubes until the 
structure of the surface inside each cube is known lfT9l . For implicit surfaces, Snyder proposed the 
globally parameterizable criterion for that purpose f28l. Plantinga and Vegter proposed the small 
normal variation condition which leads to a better meshing algorithm [24. ,25]. 

Hart et al proposed a method based on Morse theory ifTSl l23l |29l . The idea is to check when the 
topology of f{x, y, z) = a will change for a parameter a. When a changes from some initial value 
where f{x, y,z) = a has no solution to a = 0, the topology of the surface is found. Fortuna et al 
presented improved algorithms for surfaces in the projective space |[T5l[T4l . 

For a set of points on the surface, one can form the restricted Delaunay triangles and the correspond- 
ing Delaunay triangulation can be used to approximate the surface. Boissonnat and Oudot proved 
that when the sample point set satisfies certain conditions, the Delaunay triangulation has the same 
topology as the surface QIHl. Cheng et al established similar results using different strategies fV2\ . 

The CAD method proposed by Collins can be used to divide the Euclidean space into cylindrical cells 
such that the given surface has the same sign on each of the cells. Then to determine the topology of 
the surface, we need only give the adjacency information between the cells 13]|20l. Alone this line, 
new ideas are introduced to compute the topology of surfaces |[T0ll22l . 

All the above methods except the one based on CAD work for surfaces without singularities only. In 
this paper, we give a method to compute a certified meshing for implicit algebraic surfaces with sin- 
gularities. The method is a hybrid one based on the CAD approach and the marching cube approach. 
We propose a CAD based method to compute a certified meshing for the surface inside a box con- 
taining singularities and use a modified Plantinga- Vegter method to compute a certified meshing for 
the surface inside a box without singularities. Our main contribution is how to treat the singularities. 

This paper consists of three parts. The algorithms for surfaces are the main contributions. In Section 
3, a new method is proposed to compute a certified meshing for a plane algebraic curve. This section 
also provides preparations for algorithms about surfaces. There exist many methods to compute the 
topology of plane curves, e.g., |l2l|9l[T3l[l7l. Our contribution is to give an interval based method to 
compute the adjacency information and to give an ambient isotopic meshing for a curve. The method 
in 13 can also compute an ambient isotopic meshing based on root bounds of equation systems. Our 
method is based on symbolic-numerical computation, which is practically more effective. 
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In Section 4, a new method is proposed to compute an isotopic meshing for a surface. The method 
has two advantages. First, we use symbolic computation methods to guarantee the completeness 
and whenever possible use interval arithmetics to increase the efficiency. Actually, computations of 
algebraic numbers are totally avoided. The work [2', '20'] uses algebraic numbers. Second, our algo- 
rithm does not change the surface to generic positions as done in [22], which is generally expensive. 
Our method need only to project the surface once, while the algorithm proposed in ll22l need to do 
projections twice. 

In Section 5, a method is proposed to compute a certified meshing for a surface. A well-known 
technical to treat a singular point P is to find a segregating box which contains P but does not 
intersect the surface at its bottom and top faces. We f extend this concept to singular curve segments 
and give an interval based method to compute such boxes and meshes in the boxes. Another key 
ingredient is a careful analysis of the extremal points of surfaces and spatial curves. It is pointed out 
in [5 1, that the method in [22l "makes no guarantees about the geometric accuracy of the mesh, and 
it cannot be extended in a straightforward way to provide a more accurate mesh." To our knowledge, 
the method proposed in this paper is the first one to compute a certified meshing for surfaces with 
singularities. 

Algorithms in Sections 3 and 4 are implemented in Maple and nontrivial examples are used to show 
that the algorithm is quite effective for surfaces with singular points and curves. 

2 Preparations 

In this section, we give several known results and algorithms needed in this paper. Following f5l, we 
will compute a meshing with correct topology for a curve or a surface in the following sense. 

An isotopic meshing for a variety 5 C (n = 2, 3) consists of a graph/polyhedron (for n = 2, 3) 
and a continuous mapping 7 : x [0, 1] — > M" which, for any fixed t G [0, 1], is a homeomorphism 
7(-, i) from ii" to itself, and which continuously deforms W into S: 7(-, 0) = id, 7(5^, 1) = S. 

For a number e > 0, an e-meshing for S is an isotopic meshing for S, which gives an e- 
approximation for S in the following sense || P — ^{P, 1) ||< e for all P G 5^. Please note that 
isotopy is stronger than homeomorphism [5 ]. 

2.1 Real root isolation of triangular system 

A basic step of our algorithm is to isolate the real roots of a triangular system which consists of 
equations like 

= {fl{xi),f2ixi,X2),...,fnixi,X2,...,Xn)} (1) 

where /j G Q[xi, . . . , Xj] involves Xi effectively. 

We use intervals to isolate real numbers: let DQ denote the set of intervals of the form [a, b] where 
a < & G Q. The length of an interval box B„ = [ai, 61] x • • • x [a„, bn] G DQ" is defined to be 
|B„| = maxi(6i - ai). 
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In this paper, when we say a point, we mean a point P = (Ci, • • • > Cn) with real algebraic numbers 
as coordinates, which is represented by a triangular system S„ like © with P as a solution and an 
isolation box B„ for For instance ^/2 is represented by — 2 = and (1, 2). 

Now, we give a formal description of the root isolation algorithm. 

Algorithm 2.1 RootIsol(S„, B„, e). The input consists of a triangular system Tin of form ([7]), a 
box B„ G DQ", and a positive number e. The output is a set of isolation boxes for all the real roots 
of Tin = in B„ such that the length of the isolation boxes is smaller than e and any two of the 
isolation boxes are disjoints. 

A modified version of the root isolation algorithms in ifTTl l27l is used in our implementation. 

Let /(xi, . . . , Xn) G Q[xi, . . . , x„] and B„ = [ai, 6i] x • • • x [a„, 5„] G DQ". The box operation 
□ /(B„) returns an interval containing all the points {/(xi, . . . , Xn) \ ai < Xi < bi,i = I, . . . , n}. 
Furthermore, when |B„| approaches to zero, the length of interval □/(B„) also approaches to zero. 
If aj > and bi > 0, we can construct □/(B„) as follows 

□ /(B„) =/+(6i,...,&„)-r(ai,...,a„) 

where f = — f~ such that f^,f^ G Q[xi ■ ■ ■ ,x„] each has only positive coefficients and 
minimal number of monomials. For the general case, please consult [11]. It is clear that such an 
operation satisfies the following property. 

Lemma 2.2 If ^ = (^i, . . . , is not a zero of f{xi, . . . , x„) = and B„ an isolation box for 
^. Then if the length of^n small enough, the interval U f (Bn) will not contain (0, . . . , 0), which 
means that B„ has no intersections with / = 0. We denote this a^' □/(B„) ^ 0. 

2.2 Delineable polynomials 

Delineable polynomials are important in determining the topology of algebraic surfaces. Let /(xi, 
. . . , Xr-i,Xr) G M[xi, . . . , Xf] and P = (pi, . . . a point of W . We say that / has order k at 
point P, if A; > is the least non-negative integer such that some partial derivative of total order k 
does not vanish at P. And / is said to be order-invariant in a subset R of provided that the order 
of / is the same at every point of R. 

For simplification, we denote the (r — l)-tuple ( An r-variate polynomial /(x, Xr) 

over the reals is said to be delineable on a submanifold R of M*"^^ if it holds that: 

(1) the portion of the real variety of / that lies in the cylinder P x M over R consists of the union of 

the function graphs of some /c > analytic functions 6i < . . . < 9^ from R into R; and 

(2) there exist positive integers such that for every a £ R, the multiplicity of the root of /(q, x^) 

corresponding to 6i is rrii. 

Polynomial / is said to vanish identically on P if /(P, x^) = for every point P G P. In addition, 
/ is said to be degree-invariant on P if the degree of /(P, x,.) as a polynomial in Xr is the same for 
every point P G P. In this situation, the following theorem holds (see 1201 . pp. 246). 
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Theorem 2.3 (McCallum and Collins) Let f{x, Xr) be a polynomial in M[x, Xr] of positive degree 
in Xr- Let D{x) be the discriminant of f as a univariate polynomial in Xr and suppose that D{x) 
is a nonzero polynomial. Let Rbe a connected submanifold ofW~^ on which f is degree-invariant 
and does not vanish identically, and over which D is order-invariant. Then, f is delineable on R. 

The following theorem improves the above result. 

Theorem 2.4 ([8]) Let f G M{x,Xr) (r > 2) be an r-variate polynomial of positive degree in Xr 
with discriminant D{x) ^ 0. Let R be a connected submanifold of M^'~^ in which D is order- 
invariant, the leading coefficient of f w.rt. Xr is sign-invariant, and such that f vanishes identically 
at no point in R. Then, f is degree-invariant on R. 

3 Ambient isotopic meshing of plane curve 

In this section, we give an algorithm to compute an isotopic meshing for an algebraic curve. The 
main purpose of this section is to provide preliminary algorithms for later sections. We also give a 
new and fast method to compute the adjacency information based on interval arithmetics. 

3.1 Determine the topology of plane algebraic curve 

We use a graph to represent the topology of a plane curve. A topology graph is a graph W = {V,£} 
where 

• V is a. set of plane points defined by triangular systems Sj and isolation boxes Bj^: 

^ - {P^.J = P^.J), < I < S , < j < S,} (2) 

Si = {hi{x),gi{x,y)},Bij = [ai,bi\ X [ci^j^di^j] 

where ao < ai < ■ ■ ■ < as and /?j o < A,i < • • • < Pi^Si- When drawing the graph, we use 

Mij = {{ai + bi)/2, (qj + dij)/2) to represent Pij. 

• £ = {{Pi,P2)\Pi,P2 G r, such that either Pi = Pi^p,P2 = Pi+i,g or Pi = Pi^p,P2 = 
Pi.p+i}- In the first case, the edge is called non-vertical. In the second case, the edge is called 
x-vertical. We further assume that any two edges do not intersect except at the end points. 

Consider a plane algebraic curve C : g{x, y) = where g{x, y) € Q[x, y\ is a square free polynomial. 
A point Pq is an x-critical point of C if ^(Po) = 9y{Po) = 0. 

We will consider the part of C in a bounding box 

B2 = [Xi,X2]x[yi,y2]enQ^ (3) 

which is denoted as = C n B2. In the rest part of this paper, B2 is always assumed to be of this 
form. 

Let P be a point on curve C, the left (right) branch number of P, also denoted as L^{P) (R#{P)), 
is the number of curve segments of C which pass through P and are on the left (right) side of P in a 
small neighborhood of P. 
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We introduce the key concept of segregating box. A box B = [a, 6] x [c, d] G DQ^ is called segre- 
gating w.r.t. C if 

C n [a, b] x[c,c]=Cn [a, b] x [d, d] = 0. 

A curve behaves nicely in a segregating box, as illustrated by the following lemma. See Fig. I2a) for 
an illustration. 



(a) (b) (c) 

Figure 2: Curve segments inside a segregating box. 

Lemma 3.1 Let B = [a, b] x [c, d] eUQ^ be a box segregating w.r.t. C and the interior o/B contains 
no X— critical points ofC. Let C intersect the left and right boundaries ofB at points Li, i = 1, . . . ,1 
and Rj,j = 1, . . . , r respectively. Then C is delineable over R = (a, b) and the number of curve 
segments ofC inside B equals Yl\=i R-i^i^i) = Si=i ^H^i^i)- (^^^ Figure^a) for an illustration) 

Proof. Note that the leading coefficient C{x) of g{x, y) w.r.t. to y is a factor of the discriminant 
D{x) of g{x,y) as a univariate polynomial in y. Since there exist no x-critical points of C inside 
B, C{x) is not zero. Hence g{x,y) is degree invariant over R. Also D{x) = has no roots 
over R. Then, by Theorem 12.31 C is delineable over R and Cb consists of curve segments starting 
from certain Lj and ending at certain Rj. Furthermore, these curve segments do not intersect. So 
E'=i = Ei=i R#{R^)■ This proves the lemma. | 

A box B is called a segregating box for a point P on C if P is inside B, B is segregating w.r.t. C, 
and Cb \ {P} contains no x-critical points of C. See Fig. [Sja) for an illustration. It is known that 
(Theorem 5 in |12J): 

Lemma 3.2 IfB = [a, b] x [c, d] G DQ^ is a segregating box of P on C, then R^{P) and L^{P) 
are the numbers of real roots of g{b, y) = and g{a, y) = in (c, d) respectively. See Fig. ^b). 

The following algorithm computes the branch numbers. 

Algorithm 3.3 NumCur(P). V is a set of points defined by (|2]). Output R^{Pij) for < i < s — 1 
and L^{Pij) for 1 < i < s. 

1. For < i < s, if gi{x, y) has a factor of the form V{x) G Q[x\, let gi = gi/V{x). 

2. While G Qj) or G □^^([aj, djj), repeat [ai,bi] = RootIsol(/ii(j;), [ai,bi], 
{bi-ai)/2). 

3. Let i? = RootIsol((jr(6j, y), [cjj, djj], 1) and L = RootIsol((/(aj, y), [qj , djj], 1). By 
LemmaEH R#(.Pi,j) = \R\ and'L#(Pij) = \L\. (See Fig.Etb)) 

Proof of correctness. Since Bjj is an isolation box for Pij, then gi{ai, Cij)gi{ai, dij) ^ 0. By 
Lemma l2!2l the procedure in Step 2 will terminate. At the end of Step 2, giix, Cij)gi{x, dij) = 
has no real roots in [oj, hi], that is, Bj j is a segregating box for Pij. The third step is clearly true. | 
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Remark 3.4 In Step 2 of Algorithm \3. 3\ the boundary points need special consideration. If the 
boundary points are on the curve, that is, if g{cti,yi) = or (/(oj, 3^2) = 0, we will make sure that 
the following condition holds: is the only real root of gi{x, 3^i) = or gi{x, 3^2) = in [oj, bi]. 
Then, the algorithm also works. 




(a) (b) (c) 



Figure 3: Compute topology graph of a curve 

The following algorithm to compute a topology graph follows the basic idea in f7]. Our main contri- 
bution is to use interval arithmetics instead of algebraic numbers. Also, we do not need changing the 
curve to generic positions as done in |[T3llT7l 

Algorithm 3.5 TopCur(gi(x, y), B2, e). C : g{x, y) = is the curve, B2 is defined in (13), and e > 
is a number. Output a topology graph ^ = {V, £■) which is an isotopic meshing/or C-q^- Further, 
each isolation box ^ of a point in V satisfies |B| < e. 

1. Let <S = and g{x, y) = V{x)gy{x, y), where V{x) is the factor of ^(x, y) in x only. 

2. Let D{x) = Res{gv, y) be the resultant of g^ and 

3. Let V =RootIsol(S2i, B, e)URootIsol(S22, B, e), where 

H^x) = H{x)/gcd{H{x),V{x)) (4) 

521 = {Hy(x),gy{x,y)} 

522 - {Vix),g,ix,y){y-y2)iy-yi)}. 

Assume that V is of form Q- See Fig. [Ja) for an illustration. 

4. Execute Algorithm |33] to compute {1 < i < s) and R#{Pij) (0 < i < s - 1). 

5. Add an auxiliary line at 2; = {bi + aj+i)/2 and construct the non- vertical edges. See Fig|2tc) 
for an illustration. For i = 0, . . . , s — 1, execute the following steps 

(a) Let Qi =RootIsol({x — ''"^2'^^ , gv{x, y)}, B, e), where m, hi are from Arrange the 
points in Qi bottom up, we have Qi = . . . , Qi,Ui}- Set R#{Q) = L#{Q) = 1. 

(b) Let TZi be the list of points Pj ^ arranged bottom up and point ^ will be repeated 
RH^{Pi,k) times in Similarly, is the Ust of points Pi+i^t- By Lemma [XT] Ci,TZi, 
and Qi contain the same number of points. Let £j = (Li , . . . , L^^), TZi = {Ri , • • • , Rut)- 

(c) For j = 1, . . . , Ui, add {Lj,Qij), {Qij,Rj) to £. 

(d) Let P = P U Qi. Still assume that V is of form ©. 

6. Add the x-vertical edges. If ai is a root of V{x) = 0, add {Pi^k, Pi,k+i),k = 0, . . . ,Sito £. 
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7. Output the topology graph ^ = {V, £}. The isotopy map can be constructed in the usual way 

m. 



Theorem 3.6 Algorithm \3.5\ computes an isotopic meshing for C-q^. 

Proof. First, we prove that each edge e G <S represent exact one curve segment of C, and for each 
degree invariant segment of C, there exist exact one e G <5 presenting it. Hence £ and some 2/-vertical 
line decompose B2 into cylindric regions . 

It is clear that the curve C consists of two parts C„ : gv{x, y) = and V{x) = 0. The part V{x) = 
consists of straight lines x — 7^ = 0, i = 1, . . . , t, where 7j are the real roots of V{x) = 0. To 
determine the topology of C, we need only to find the topology graph '^^ of and then to add the 
Unes X — 7i = to ^^t,. So we may consider only. 

From Steps 2-4, we know that V contains all the x-critical points of the curve C„ and the boundary 
points which are the intersection points of C and the boundaries of B2. In Steps 6 and 7, we add 
auxiliary points Qi to V. Since points in Qi are not critical points of C, we have R^{Q) = L#{Q) = 
1 for Q G Qi. This makes sure that all the edges {Lj, Qij) and {Qij,Rj) are distinct. 

Let Bi = (oj, Oj+i) X [3^1, 3^2], i = 0, . . . , s — 1. We need only to show that and have the 
same topology in Bi. Let Si be the interval (ai,aj+i). Then D{x) does not vanish on any point 
of Si. As a consequence, gy{x,y) must be degree invariant on Si. By Theorem 12.31 gv{x,y) is 
delineable over Si and '^^ is obtained by replacing a curve segment of in Bi by a line segment 
with same end points. It is clear that Cy and '^^ have the same topology. We are going to make 
explicit the isotopy from £ to Cb2- Let = {V,£) be a topology graph for curve Cb2> and V of 
form ©. Let Pjj = {ai,pij) be of form ©. Let Qij = {n^pij) where n = ^^,Pi,j = "''^^^'-^ 
Then ^ decomposes B2 into cylindrical regions UijRi^k, where Rij is bounded by [rj, Tj+i] in the 
x-direction and by /i = {Qi^u, Qi+i,v) and /2 = (Qi.si Qi+i,t) for ceratin u, v, s, t. Note that i?j ^ 
could be a triangle or a quadrilateral. 

First, we consider one cylindrical region i?^ ^ defined as above. Let ei = {Pi,u, Pi+i,t,)) and 62 = 
(Pi^s, Pi+i,t)- Without loss of generality, assume ai < 02, < Pi,2- According to the correctness 
prove of Algorithm 13.51 gv{x, y) is delineable over [ai, a2]> we can find two root functions 6i{x) of 
g^ on [01,02] corresponding to the two curve segments C(ei) and C(e2). Denote y = 5i{x),x G 
[01,02] to be the definition functions of line segments Cj and y = ipi{x),x G [Ti,r2] to be the 
definition functions of line segments /j. Consider the maps: 



Fi : ([01,02] X M) X [0,1] ^ 



defined by 



{x,XSi{x) + {l- XS2{x)),t) 
{x,X{t9i{x) + {l~t)Siix)) + {l 



X){t92{x) + {1 - t)S2ix))) 



and 



F2 ■■ ([ri,T2] xM) X [0,1] 



[01,02] X M 



defined by 



{x,Xipi{x) + (1 - X^p2{x)),t) 
ix',Xitdiix') + il-t)^i{x)) + il 



X){t62ix') + {l-t)ip2ix))), 
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where x' = ai + 



X — T\ 

T2-ri 



(a2 - ai)- 



The map F2 is a homeomorphism from [ri , x M to [ai , Q!2] x M and Fi is a homeomorphism from 
[ai, 02] X M to itself. So the composed map Fj := Fi o F2 is a homeomorphism from [ti, r2] x M 
to [ai, 02] X M and it deforms fi to C(ej) continuously. Extend this map to x [0, 1] by setting 
it to be the identity map outside we obtain an isotopy from line segments /i U /2 to the curve 
segments C(ei) U C(e2). 

Now we consider the whole topology graph For each cylindrical region Ri^k, we can construct an 
isotopy Fj jfc as above. Consider the following map: 



Note that Fjj|R, ^.nH„,„ = -^«,t;|i?,,,nR„,„, and Fij|^^^^.n(iR2\B2) = for all i,j,u,v. (?^,F) is an 



As a consequence of the above proof, we have 

Corollary 3.7 Let G = (V,£) be a topology graph of the curve C-q^ obtained by Algorithm \3.5\ 
Then all the singular points of 'Sire in V and g{x,y) is y -degree invariant over the intervals 
{ai,ai+i),i = 0, . . . ,s - 1. 

When computing the topology of a surface, we need to introduce the concept of extended topology 
graph. An extended topology graph associate with a box B2 is a triplet £Q = {£V, ££, £C} where 



{£V,££} is a topology graph and £C = {{P^, P2, Ps)] B, G £V, (Fi,F2), (P2,i^3), {P3,Pi) 



£ ££} is a set triangular cells in B2. We further assume that the cells in £C are disjoint except on 
their edges and provide a cover for B2. 

We can obtain an extended topology graph of a curve from a topology graph by adding more auxiliary 
points and edges. 

Algorithm 3.8 ETopCur((7(x, y), B2, p) The input is the same as Algorithm \3.5\ The output is an 
extended topology graph o/Cb2- (See Fig. \3\c)for an illustration) 

1. Let?^ = {V,£] = TopCur(5(2;,y),B2,p). 

2. Let £V = V. For z = 0, . . . , s, add (oj, y\) and Fj = (oj, 3^2) to £V if they are not in it. 

3. For points Pij,j = 0, . . . , Sj, let [oj, hi] x [cij,dij] be the isolation box for Pij. Add Nij = 
[ai, (dij + Cjj+i)/2), j = 0, . . . , Sj — 1 to £V. We still assume that £V is of form (|2]). 

4. Let ££ = £. For j = 0, . . . , sq — I, add the edges [Pqj, Nqj), (Nqj, Pqj^i) to ££. 

5. For each < i < s - 1, add the edges (Fj^o, -Pj+1,0) and (Fj,^., Fj+i^^.^J to ££. Then 
the edges in ££ divide the rectangular region B = [ai,ai+i] x [3^1,3^2] into triangular and 
quadrilateral regions. We will subdivide these regions into triangular regions such that each 
point in £V is the vertex of at least one triangles. 



F : X [0, 1] 



denoted by 




(5) 



isotopy for ■ 



9 



6. Let EC = 0. For any two adjacent edges ei = (Pi,i, i^2,i), 62 = (-Pi,2, -^2,2) inside B, execute 
Steps 7 and 8. 

7. If Pi^i 7^ Pi^2> there exists one point A'^i added before in Step 3 between Pi^i and Pi^2- 
Furthermore, 

• If ^2,1 / ^^2,2. there exists one point between P2,i and ^2,2- ^1,1, ^1,2, -^2,2, ^2,1 
form a quadrilateral region. We can divide the quadrilateral region into four triangles. 
Add the edges (^2,1,^^2), (A^2, ^2,2), (^1,2, N2), {Ni,N2), (iVi, ^2,1) to ££. Add the 
triangles iVi, P2,i), (A^i, ^2,1, iV2), (iVi, A,2, iV2), (Pi,2, iV2, ^2,2) to ^C. 

• If ^2,1 = ^2,2 = P, Pi,i) Pi.2,P form a triangular region. We can divide the triangular 
region into two triangles. Add the edges {Pi^i,Ni), {Ni,Pi^2), {^i,P) to ££. Add the 
triangles (Pi,i,iVi,P),(iVi,Pi,2,P) to ^:C. ' 

8. If Pi^i = Pi^2 = P> then there must exist a point A'^2 added before in Step 3 between 
P2^i and P2.2- P, P2,25 P2.1 form a triangular region. We can divide the quadrilateral region 
into two triangles. Add the edges (P2,i,iV2), (A'^2,P2,2), (P, ^^^2) to ££. Add the triangles 

(P,P2,l,iV2),(P,P2,2,iV2)t0fC. 

9. Output £g = {£V, ££, £C}. 

Remark. The purpose to add points Nij in Step 3 is to make sure that topology representation for 
surfaces possible. These points has similar function as the the auxiliary points added in Step 6 of 
Algorithm l3.5l Hence, they are also called auxiliary points. Figure[3]is an extend topology graph of 
the curve G{x, y) = x ■ y ■ (IGx^ + 16y^ - 49) = 0. 

Let £Q = {£V,££,£C} be an extend topology graph of and e = (Pi,P2) G ££. If e cor- 
responds to a curve segment of Cb2, we use C(e) to represent the corresponding curve segment; 
otherwise, we use C(e) to represent the line segment P1P2. Let /(e) = C(e) \ {Pi,P2}. For a 
c G £C, we use R{c) and /(c) to denote the cell and interior of the cell represented by c respectively. 



3.2 Compute e-meshing for plane curve 

The meshing given in Section [3?T] has no guarantee of precision. In this section, we will show how 
to compute a meshing for a curve to any given precision. 

Let = {V, £) be a topology graph for a curve C inside a box B2 defined in (O. Assume that V is 
of form Q. Consider the two disjoint regions S2 and N2 of B2: 

S2 = UiS'2,N2 = UjN^, where 

Si = (ai,60 X [3;i,3^2],i = 0,...,s (6) 
N^2 = h,aj+i] X [3^1,3^2], J = 0, . . . , s - 1. 

Then, C S2 U N2 and is smooth in N2. 

The idea of our algorithm is to determine the topology of the curve in the region S2 with Algorithm 
13.51 to determine the topology of the curve in the region N2 with a modified marching cube method of 
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Pantinga-Vegter II24II . and to compute the adjacency information on the border lines 
We could use the marching cube method in N2 because C has no singular point in it. 




(b) 



Figure 4: Nice boxes: (a), (b). Boxes in (c), (d) not nice, c at u- 

^ v/>vy v/>vy Figure 5: Meshmg curve segments 



In order for the above idea to work, we need to modify the Pantinga-Vegter method such that each 
output box contains only one curve segment of C, as shown in Fig. Ill^a) and (b). Such boxes are 
called nice boxes. 

The original Pantinga-Vegter method could output a box containing two curve segments and this will 
cause problems when the box is near a singular point, as shown in Fig. HJc). A point is called a 
y-extremal point of curve C if C achieves a local extremum value at this point in the y-direction. 
Pantinga-Vegter's method could output a box shown in Fig. I5d). 

To make the process precise, we introduce the following definition. An e-meshing graph of a curve 
C is a triplet Ai = {V, £, B} where {V, £) is a graph whose vertices are with rational numbers as 
coordinates and whose edges are the meshes for C; ;B is a set of nice boxes and segregating boxes 
of singular points of C such that for each e G <S, there exists a Be € ;B with the property: |Be| < e 
and C n Be is a connected curve segment of C (See Fig. [S]). In Fig. 12 a), V = {Ni,N2},e = 
[Ni, N2),'Be = ABCD forms a meshing graph for curve segment C(e) = P1HP2. 

It is easy to show that an e-meshing graph for a curve C provides an e-meshing for C according to 
the definition given in Section 2. 

Algorithm 3.9 MPV2((7(x, y), B2, e). Input: C : g{x,y) = is a curve with no x-critical points 
and no y-extremal points in box B2. Output an e-meshing graph ^ = {P, E, B}for Cb2- 

We need only add some extra criterions for the boxes: (1) For each edge {A, B) of B, if G 
Ug{{A, B)) and g{A)g{B) > 0, we continue to subdivide B. (2) For each box B, if |B| > e, 
we continue to subdivide B. 

Since C has no x-critical points and no y-extremal points in box B2, a box like the one in Fig. |4td) 
does not exist and the algorithm will terminate. 

Now, we can give the meshing algorithm for curves. 

Algorithm 3.10 ATopCur(5'(x, y), B2, e). The input is the same as that of Alsorithm \33\ Output 
an e-meshing graph for g{x,y) = 0. 

1. Execute the first four steps of Algorithm 13.51 with input (g{x, y), B2, e). We need to modify 
Algorithm 13. 5l as follows: Let gu{x, y) = gv{x, y)/U (y) where U{y) is the gcd of the coeffi- 
cients of gv{x, y) as a univariate polynomial in x; and use H{x) ■ Res((7u, ^^,y) as the new 
H{x) in @. 

We need V{x), gi,{x, y), and '^i = {Vi,Bi} from Algorithm [331 where Bi is the segregating 
boxes for the points in Vi. 
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2. Compute ^2 = {^2, £2, S2}=MP\2{g{x, y) ,^2, e). The modification in Step 1 makes sure 
that C has no x-critical points and y-extremal points in N2. 

3. Compute the connection between the boxes computed by Step S 1 and Step S2. (Fig. |5tb) 
shows how to mesh C near a singular point O with B = ABCD as its segregating box. Fig. 
Ob) provides a global picture for meshing a curve.)For i = 0, . . . ,s, consider the adjacency 
information on the border lines x = ai,bi. We only consider x = hi. For each P ^ V and its 
segregating box B = [a, b] x [c, d] G B, do the following 

(a) Let Efc = [b, c^] x [e^, fk] G B2 be the boxes satisfying BnEj. / and g{b, ek)g{b, fk) < 
0, where ek = minjefc, c} and fk = maxj/^t, d}. As a consequence, C passes through 
these Ejt through the interval [6,, 6j] x [sk, fk]- 

(b) Let Q = ((a + 6)/2, (c + d)/2), nik = (4 + A)/2. 

(c) Add the edge e = {Q, (6, 771^)) to M. Add Be = B to A^. 

4. Add the meshes for the straight lines defined by V{x) = 0. 

5. Output the meshing graph A^. 




Figure 6: e-meshing for a curve 



Theorem 3.11 Aleorithm \3. lOU erminates and computes an e-ambient meshing for Cb2- 

Proof. By Lemma [331 H S2 is the isotopic meshing of C n S2. Marching cube method compute 
the isotopic meshing of C n N2. Hence, is a isotopic meshing for Cb2- Furthermore, for each 
e G <5 n N2, C(e) C Be, I Be I < e and each part of around the singular point P is contained in 
the segregating boxes Bp, |Bp| < e of P. Therefore, for any point P ^ M, F{P, 1) and P are in 
the same box Bg with |Be| < e, so || F{P, \) — P ||< e. This gives aproof of Theorem |3.11| | 

In order to compute the e-meshing for surfaces, we need to add more information to the e-meshing 
graph. Let M. = {V, £, B} be a meshing graph for a curve C. Then an extended meshing graph 
£Ai = {£V,££,£C} for C in B can be defined similarly as the extended topology graph. The 
difference is that £C provides a triangular decomposition for B. 

The following algorithm computes an extended meshing graph. 

Algorithm 3.12 METopCur((7(x, y), ^1, ?^2)- '^1 = {'PijBi} where Vi is the set of points on 
C : g{x,y) = of form ([2]) and Bi is the set of their segregating boxes. '^2 = {'P2,£2,B2} is an 
e-meshing of the curve C in N2 defined in Output an extended meshing graph for C in Bi U B2- 
(Fig. ^c) is the extended meshing graph for the box in the center of Fig. ^a) and its surrounding 
boxes.) 

SI Let£V = $,££ = $,£€ = 0. 
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52 For any B = [a, b] x [c, d] & Bi, it is the segregating box of one point P e Vi. Compute the 

extended topology in Bi. 

1. Compute 

{Qi,...,Qi} = RootIsol({j; - o,5(a,y)}, [c,(i],l) 
{Ti,...,Tr} = RootIsol{{x-b,g{b,y)},[c,d],l)- 

Denote Ai,i = 1, . . . , 4 to be the four vertices of B. Denote Bj,j = 1, ... ,p to be 
the points on the edge of B which are the vertices of boxes adjacent to B.(Note that if 
/ > l(r > 1), there exists some point between Qi and Qi^i{Ti and Tj+i)). 

2. TV = %. Add points Qi,Tj,Ak,Bp into TV. 

3. Denote the sets of points L = {Li, . . . , s} and R = {Ri , . . . ,Rt} where L is the points 
in TV which are on the left edge of B sorted from bottom to up and R the points in TV 
which are on on the right edge of B sorted from bottom to up. 

4. Add point P and all points in TV into £V. Add edges (P, Li), {P, Rj) into ££. Add 
triangular cell (P, Lj, Lj+i), (P, P^, P^+i) and (P, Li, Pi), (P, L^, Pt) and (Li,Li+i), 
{Rj,Rj+i) to SC. 

53 For any boxes B G ^62. Compute the extended topology in B. For any line segment e £ £2 with 

B = [a, 5] X [c, d] G B2 containing it. Doing the following operations(There are six conditions 
that e divide B into two parts, see fig ID We can distinguish them according to V2 and B2. 
Here we consider the condition (a), the other conditions are dealt with in the similar way). 

1. Compute Q = RootIsol({x— a,5(x,y)},B,e/4) and T = RootIsol({3; — 6, y)}, B, 
e/4). Obviously, Q and T both contain only one point. We still call them Q and T. 
Denote Ai,i = 1, . . . , 4 to be the four vertices of B. Denote Bj,j = 1, ... ,p to be 
points on one edge of B which are the vertices of boxes adjacent to B. 

2. TV = 0. Add points Q,T,Ai, Bj into TV. 

3. Add all points in TV into £V. Similar to the forth step in Step S2, we can decompose 
B into triangular cells and insert these cells into £C, and insert corresponding edges into 
££ such that each point in TV connects to at least another point in this set. 

54 Output fTW = {£V,££,£C}. 

4 Topology of surface 

In this section, an algorithm will be given to compute a polyhedron with triangular faces, which is 
isotopic to a given surface. 
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4.1 Outline of the algorithm 



We use a polyhedron with triangular faces to represent the topology of a surface. A topology poly- 
hedron is a triplet 3^ = {SV, S£, ST} where SV, S£, and ST are defined below. 

• SV is a set of 3D points determined by a triangular system and an isolation boxes ^ijy. 

SV = {P,j,fc, < z < s, < J < s,, < fc < <,j} 
Si = {hi{x),gi{x,y),fi{x,y,z)} (7) 

where P^j- ^ = {ai, Pij,jij^k) satisfy ao < • • • < "s, A,o < • • • < and 7^^- o < • • • < 
HJA J - Point is said to be lifted from the plane point Pjj = {ai, f3ij). Pij is said to be 
the projection of Pij.k- 

• S£ = {{Pi,P2)\Pi,P2 G SV, such that either Pi = P,,„,„,P2 = P^+l,p,q or Pi = 
Pi,u,v,,P2 = Pi,u+i,t}- We further assume that any two edges do not intersect except at the 
end points. 

• ST = { (Pi , P2 , P3 ) I Pi , P2 , P3 G SV} such that its three edges are in S£. We further assume 
that any two faces do not intersect except at the edges. 

Let S : f{x,y,z) = be an algebraic surface, where f{x,y,z) G Q[x,y,z] is square free. A 
point Pq is a critical point of S if /(Pq) = fz{Po) = 0. Write / as a univariate polynomial in z: 

f{x, y, z) = fd{x, y)z'^ H h fo{x, y). fd{x, y) is called the leading coefficient of f{x, y, z). We 

further assume that 

fd{x,y) = --- = fo{x,y) = have no common zeros. (8) 

Geometrically, this means that S does not contain a line parallel to the z-axis. We will consider 
surfaces that do not satisfy this condition in Section 14771 

Similar to the case of algebraic curves, we will consider the topology of 5 in a bounding box 

B3 = [Xi,X2] X [3^1, 3^2] X [Zi,Z2] G DQ^. (9) 



Let 

df 

D{x,y) = Res(/,/,z) (10) 
oz 

G{x,y) = sqifme{D{x,y)f{x,y,Zi)f{x,y,Z2)) (11) 

where sqrfree(P(x, y)) is the square free part of P{x, y). The plane curve G{x, y) = is called the 
projection curve of S. 

To determine the topology of a surface is to find a topology polyhedron with the same topology as 
the surface. We first give an outline of the algorithm, which consists of four main steps. 

SI Compute an extended topology graph £Q = {£V,££,£C} of the projection curve of S in B2 
defined in ([3]). 
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52 Determine SV. For any P G £V, determine the intersection points of S and the line segment 

53 Determine S£. For each edge e G ££, compute the intersection of S and the cylindrical surface 

patch /(e) X [Zi, Z2], which are delineable curve segments of S whose end points ai^e in SV. 
We will use line segments in S£ to represent these curve segments. See Fig. |7] 

54 Determine SF. For each c G £C, compute the intersection of S and the prism /(c) x [^1,^2], 

which are delineable surface patches of S whose edges are in S£. We will use triangular faces 
in ST to represent these surface patches. See Fig. [8] 

4.2 Theoretical preparations for the algorithm 

In the outline of the algorithm given in the preceding section, Step SI has been solved in Section [3?T] 
Step S2 can be solved with Algorithm Rootlsol. We will explain Steps S3 and S4 below. 

Roughly speaking, Step S3 is to determine the topology of the spatial curve defined by f{x, y, z) = 
G{x, y) = 0. The following result, which is a consequence of Theorem 12.31 allows us to determine 
the singularities of this curve easily. 

Lemma 4.1 Use the notations introduced above. For each edge e = {Pi, P2) £ ££, f{x, y, z) = 
is delineable over /(e) = C(e) \ {Pi, P2}. 

Proof. Let C : G{x, y) = be the projection curve of S and D{x, y) the discriminant of / w.r.t. z. 
Since /(e) is a continous curve segment of C, G is order-invariant on /(e). From (ITOl ). D{x, y) is 
order-invariant on /(e). Since condition dSJl holds, / does not vanish identically on any point of 
xy-plane. So, / does not vanish identically on /(e). Now, we will prove that / is degree-invariant 
on /(e). It is clear that all the singular points of C are in £V. Then /^(x, y) is either identically zero 
on /(e) or does not vanish on any point on /(e). So we can conclude that fd{x, y) is sign-invariant 
on /(e). By Theorem 12.41 / is degree-invariant on /(e) . By Theorem 12. 31 / is delineable on /(e). | 

As a corollary, we have 

Corollary 4.2 For e = (Pi, P2) G ££, the intersection ofS and /(e) x [^1, ^2] consists of disjoint 
curve segments of S whose end points are in SV. 

These curve segments together with their endpoints are called the spatial cylindrical curve seg- 
ments (SCCS) of S lifted from e. 

To determine the edges of the topology polyhedron, an SCCS with end points Pi and P2 is repre- 
sented by the line segment e = (Pi, P2). S£ is the set of these line segments. For an edge E G S£, 
we use S{E) to denote the corresponding SCCS of S. 

Let Pij^k G SV and e = {Pij,Pu^v) £ ££■ We use i^{Pij^k, e) to represent the number of SCCSes 
which have Pij^k as an end point and are lifted from C(e). We use #(e) to denote the number of 
SCCSes lifted from e. Define i^{Pu,v,w, e) similarly. As a direct consequence of Lemma l4~n the 
following equation 

#(^) = Y.*(^^'^'k,e) = J2#{Pu,v,y.,e) (12) 

k w 
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holds for each e = {Pij, Pu,v) G ££■ (See Fig. IT]) 

In Step S4, we find the surface patches lifted from a triangular cell c £ £C by identifying their 
boundaries which are SCCSes of S. As a consequence of Theorem 12. 31 we have 

Lemma 4.3 Let c £ SC. Then f{x, y,z) =0 is delineable over S = I{c). 

Proof. For any P = (a, /3) G S, / is degree-invariant and does not vanish. The discriminant D{x,y) 
of / does not vanish on P. So D is order-invariant over S. By Theorem 12. 31 the lemma holds. | 

Lemma 4.4 S n (/(c) x [Zi,Z2]) consists of disjoint surface patches whose edges are SCCSes and 
whose vertices are points in SV. These surface patches with their edges and vertices are called 
triangular surface patches (TSP) lifted from c. 

Proof. By Lemma 1431 the intersection of S and /(c) x [Zi,Z2] consists of disjoint surface patches. 
The edges of a surface patch s are the intersection of S and /(cj) x [Zi,Z2],i = 1,2,3, where 
Ci are the three sides of c. As a consequence, the edges of these surface patches are SCCSes. If 
c = (Pi,P2,/'3), the vertices of an intersection surface patch are the intersection points of S and 
Pi X [Zi,Z2],i = 1, 2, 3. As a consequence, the vertices Qi, Q2, Q3 of a triangular surface patch 
are points in SV lifted from Pi , P2 , P3 respectively. I 

It is clear that the TSPs are the intersection of C(e) x [Zi,Z2] and S. 

For a cell c = {Pij, Pu^v, Ps,t) G £C and an edge E = {Pij^k, Pu,v,w) £ S£ Ufted from the side 
e = {Pij,Pu^v) of c, we use #(c) to denote the branch number of TSPs lifted from P(c) and 
#(£^, c) to denote the number of TSPs which pass through S{E) and lifted from P(c). Notations 

v,w , Ps,t 

,i), c) and #((Pij,fc, Ps,t,i), c) can be similarly defined. As a consequence of Lemma 
1441 for c = {Pij,Pu,v,Ps,t) e £C, we have 

= E ^) = E #(^2, c) = ^ #(P3, c), (13) 

E-i 

wherePi = (Pij, fci, /'u,!,,^)' -^2 = {Pu,v,k2^ Ps,t,kz)^ = (-Pi,j,fci> -P^.t.fcs) for ^11 possible ki,k2,h. 
(See Fig. [Hi 




Figure 7: Mesh SCCSes Figure 8: Mesh TSPs 

4.3 The algorithm 

Following the analysis in the preceding section, we now give the algorithm to construct a topology 
polyhedron for a given surface. 
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Algorithm 4.5 TopSur(/(3;, y, z), B3). S : f{x,y,z) = is the surface satisfying condition iHSljf 
and f is square free. B3 is defined in dP]). Output an isotopic topology polyhedron 0^ for S-q^. 

1. Compute the projection curve C : y) = as in (fTTI) . 

2. Compute the extended topology graph: EQ = {£V, ££, £C} of with Algorithm I3.8[ 
where B2 is defined in 

3. Compute SV. For any Pi^j e £V, use Algorithm 14.71 with input (/, B3,Pij, 1) to compute 

Pi,j,k- 

4. Compute S£. Let 5f = 0. 

(a) For each t G £V and e G <5<? with P^,* as an endpoint, use Algorithm |4.9| to compute 

(b) For any e = (Pjj , Pu,v) G , let Li = {Pijfl, . . ., Pi,j,s, , ) such that point Pij^k repeats 
#{Pi,j,k, e) times. Similarly, define L2 = {Pu,v,o, • • • , ^'«,i;,s„,J- 

(c) By (Ell, |Li| = IL2I = m. Let Li = (Pi,...,P^) and L2 = {Qi,...,Qm)- Add 
(Pj, Qi) to 5<S. See Fig. |7]for an illustration. 

5. Compute 5 J^. Let cSJF = 0. 

(a) For each cell c G and E G 5<S lifted from a side of c, compute c) with Algo- 
rithm gJI] 

(b) Let ei, 62, 63 G be the three sides of c. Let 5, be the sequence of edges in S£ lifted 
from Cj ordered bottom up and an E is repeated c) times in the sequence. 

(c) By (ini), |5i| = 1521 = IS-sl = t. Let Si = {Ei^k, /c = 1, . . . Then the three line 
segments Ei^k, E2,k, Es^k should form a triangle / = (Pi_fc, P2,k,Pz,k)- Add / to SF. 
See Fig. [8]for an illustration. 

6. Output = {SV, S£, SJ^}. The isotopic map can be computed as usual ||25]| . 
Theorem 4.6 Algorithm \4. 51 computes an isotopic meshing for S-q^. 

Proof. First, we prove the algorithm compute the correct topology of given surface. Note that with 
the auxiliary points added in Step 6(a) of Algorithm 13.51 and Step 3 of Algorithm 13.81 the edges in 
Step 4(c) and the faces in Step 5(d) are mutually different. Thus, we have a well-defined polyhedron. 

The extended topology graph EG divides the rectangle B2 into triangular cells. We need only to 
show that for each edge e G ££ and each cell c G <?C, ^ and S have the same topology on 
C(e) X [^1,^2] and C(c) x [^1,^2] respectively. 

For e G ££, from Step 4 the SCCSes of S on the cylindrical surface Si = (7(e) x [^1,^2] do not 
intersect except at the end points. By Corollaries I4.2l and ([T2l) . the edges of S£ are the line segments 
with the same end points as those SCCSes. Then, the plane graph ^ on e x [Zi , Z2] and S on 5i 
have the same topology. See Figure|7]for an illustration. The spatial curve segments are presented by 
line segments. With similar arguments, we could show that the part of on c x [^1,^2] and S on 
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C(c) X [Zi, Z2] have the same topology. See Figure [8] for an illustration. This proves the topology 
correctness of the algorithm. 

Then we prove the topology polyhedron is a isotopic meshing of the given surface. 

The extended topology graph £G = {£V, ££, £C} for the curve decompose B2 into triangular 
cells. According to Theorem 13.61 £Q and C are isotopic and we can construct a homeomorphism F 
from to itself that deforms ££ to C continuously: 

F : X [0, 1] R^. 

Let ^ = {SV,S£,ST) be a topology polyhedron for a surface Sbs which decomposes B3 into 
cylindrical regions in a similar way as described in the proof of Theorem l3.6l Extend F to x [0, 1] : 

Ti = {Fix,y),z) -.W" X [0,1] ^ Ml 

The inverse transformation T^^ of Ti deforms all SCCSs of S into planes {C(e) x M, e G ££} 
which are perpendicular to the xy-pane. Denote Si to be the surface T^^{S). We need only to prove 
that Si and ^ are isotopic. 

We can construct a homeomorphism T2 from to itself similar to that give in the proof of Theorem 
l3.6l to deform the z direction such that T2{ST, 0) = ST and T2{ST, 1) = Si. 

The transformation T = Ti o T2 is a homeomorphism from R'^ to itself which deforms ST to S 
continuously. | 

We implemented Algorithm I4.5l in Maple. Two groups of experiments are done for the following five 
surfaces with singularities. 

51 : /i = a;-* + / + 2* - - y2 - - x'^y^ - x^z'^ - y^z'^ + 1 = 0, B3 = [[-1.5, 1.5], [-1.25, 1.25], [-2, 2]]. 

52 : /2 = -1 + {27/2)z'^y^x^ - {27/2)x^y'^ - Qx^z^ - {27/2)y'^z^ + 3a;^ + 3z^ + (27/4)y^ - 3a;'' - (243/16)2/" - 
Zz'^ + x^ + (729/64)^" + z'^ + (27/4)x*?/2 + Zx^^z"^ + {2AZ/l&)x^y^ + ix^z-^ + (243/16)2^1/'' + (27/4)2*^^ -x^z^- 
(9/80)2/^2^ = 0, B3 = [[-2, 2], [-2, 2], [-4, 4]]. 

S-i : f3 ^ -2^" + 2y^z^ + y^ + z'^ - 2z^ + x^ + Sx'^y^ - 3s* + Sx^y'^ - 6x^y^ + 3x^ + y** = 0, B3 = 
[[-2, 2], [-2, 2], [-2, 2]]. 

54 : /4 = x^y" + y^z^ + z^x^ - 7xyz/2 = 0, B3 = [[-2, 2], [-2, 2], [-2, 2]]. 

55 : /s = 16-22:^22-82^+4x^-2:^ + (l/4)a;*^+a;'' + y'' + ?/^2:^ + 2'' + 2^2;^-22:^y^ + 2?/^2^-8a;2-8y^ = 0, 
B3 = [[-2,2],[-3,3],[-6,6]]. 

The first experiment is to compute an isotopic polyhedron for the surfaces without considering preci- 
sion. The timings are given in the second row of Table [U Two of the polyhedrons are shown in Fig. 
|9] In the second experiment, we continue to subdivide the intervals between [Xi, X2] to compute a 
more accurate meshing. The results are given in Fig. 1. The timings are given in the third row of 
Table [T] 7^ Mesh in the fourth row gives the number of meshes in these meshings. Considering that 
implementations in Maple are generally slow due to overhead costs, our algorithm is quite effective. 
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Figure 9: Topology polyhedrons for surfaces Si and Figure 10: Isolation intervals 



TYPE 


Si 


S2 


S3 


Si 


55 


Topology 


0.544 


0.816 


0.760 


0.684 


1.280 


Meshing 


11.7 


11.8 


22.0 


51.1 


92.0 


#Mesh 


1472 


1612 


3032 


3658 


5456 



Table 1: Timings on a PC with Linux OS, 3.00G Core 2Duo CPU, and 2G RAM. 

4.4 Segregating box for a point on S 

Assume that SV is of form Then Bj j in Q is an isolation box for Pij and Bj j is an isolation 
box for Pjj,fc. It is clear that 

f{ai,Pij,eij^k)f{ai,Pi,j,dij^k) (14) 

The isolating box ^ij^k of Pi.j,k is called a segregating box if f{x, y, z) does not intersect with the 
top and bottom faces of ^ij,k- Due to (fT4)) . when sufficiently subdividing Bj Bj ^ will become a 
segregating box. This leads to the following algorithm. 

Algorithm 4.7 SegBoxP3(/(x, y, z), B3, P, e) where S: f{x, y,z) = is the surface, B3 defined 
in dP]), P a plane point defined by S2 = {h{x),g{x, y)} and an isolation box B, and e > 0. Output 
the set of points {Pi} on S lifted from P, segregating boxes for Pi, and a new segregating box B of 
P. 

1. Let{Bi, . . . ,BJ = RootIsol(S3,Bx [Zi,Z2],e), where S3 = {h{x), g{x,y), f{x,y, z)}. 

2. Let Bj = B X [ej, /»] be the isolation box for Pi on S. 

3. Let T] = e. While G □/(B x [cj, Cj]) or G □/(B x [/j, /j]) for some /c G {1, . . . , s}, repeat 

T] = rj/2 and B := RootIsol(i;2, B, rj). 

4. Output the points Pi defined by S3 and Bj, and the new B. 

In Step 3, if /(oj, Aj , 3^i) = or f{ai,Pij,y2) = 0, then we need to use the minimal circle method 
introduced in lITOl to find a segregating box in order for Lemma l4T8] to be true at this boundary point. 

4.5 Compute number of SCCSes adjacent to a point 

Let Pij^k be a point lifted from point Pi j and e = {Pij,Pu,v) G ££■ We will show how to compute 
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For any point P on the projection curve C : G{x, y) = and a segregating box B = [a, b] x [c, d] of 
P, C intersects only with the vertical boundaries of B. 

For an edge e £ £, consider the right boundaries of B. We denote the intersection point of C(e) with 
line X = b as Qe and 

[b,b] X [Ue,Ve] (15) 

is an isolation interval for Qg on line x = b, which is called the isolation interval of C{e). See Fig. 

m 

Lemma 4.8 Use the above notations. If Bj j ^ is a segregating box for Pij^k <^f^d 'S is deline- 
able over /(e), then i^{Pij^k, e) equals to the number of solutions of the triangular system S/j = 
{G{bi,y), f{bi, y, z)} in the interval box [tig, Ve\ x [cjj-^fc, /i,j,fc]- Geometrically, this is the number of 
intersection points of the line segment {x = bi,y = 'ji, Cjj-^fc < z < fij^k} '^^d the surface S where 
{bi,^i) is a point on G{bi,y) = 0. See Fig. \Tl\for an illustration. 

Proof. From Algorithm 13.31 each SCCS passing through Pij^k and projecting to C(e) must pass 
through the rectangle [6;, bi] x [ug, Ve] x [Zi, Z2]. Since Bjj ^ is a segregating box, these SCCSes 
must intersect with the the rectangle TZ = [bi, bi] x [ue, fg] x [&i,j,k, fi,j.k\- Further, each SCCS can 
intersect with the rectangle only once since these SCCSes are delineable by Lemma l4~n Note that 
the number of solutions of the triangular system S/j is the number of intersections of the SCCSes 
and the rectangle TZ. I 

Remark. Similarly, we can compute the number of the SCCSes on the left hand side of the point 
Pi,j,k by computing the number of solutions for {G{ai,y) = 0, /(a^, y, z) = 0}. When G{x, y) = 
contains vertical lines, we can compute the number of SCCSes passing through Pij^k and projecting 
to these lines by solving {G{x, w),f{x, w, z)} for w = Cij and w = dij respectively. 




Figure 11: Compute #{Pij^k^ e) Figure 12: Compute #(5", c) 

We now give the following algorithm to compute the number of curve branches. 
The following algorithm is based on Lemma [481 

Algorithm 4.9 NumSCCS(/(3;, y, z), Pij.k^ e) S : f{x, y,z) = is a surface delineable over /(e), 
Pi,j,k £ 'SV is of form (O, and e G ££ is an edge with Pij as an end point, where Pij is the 
projection point of Pij^k- The output is #{Pij^k: e). 

1. If e is an j;-vertical line segment above Pij in the y-direction, then form the triangular sys- 
tem 1:22 = {gi{x, dij), f{x, dij,z)} and let Q = RootIsol(S22, [oi, bi] x [eij^k, fi,j,k]' !)• 
Output#(Pij-fc,e) = \ Q\. 
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2. If e is not an x-vertical line segment, we need to compute tiie isolation intervals defined in 
([Tsl l. We only consider the rigiit branciies. Let TZ = RootIsol((/j(6j, y), [cij,dij], 1) wiiere 
r = R#{Pij) = \Tl\. By Lemma ll!2l one of intervals in TZ is the isolation interval [ue, Vg] for 
e. 

3. Let S21 = {gi{bi,y), f{bi,y, z)} be a triangular system in y and z and Q = RootIsol(S2i, 
[ue,Ve\ X [eij,fc,/ij,fc],l)- Output # (Pi j,fc,e) = \ Q\. See Fig. [m for an illustration. 

If there exist no SCCSes originating from a point, it is an isolated singularity. 
4.6 Compute number of TSPs adjacent to an SCCS 

We compute the number of TSPs originating from an ii^ G S£. That is, for an ii^ = {Pij,k, Pu,v,w) £ 
S£ and a c G <SC with e = {Pij,Pu,v) as an edge, we will compute #(£', c). 

Use the notations in Algorithm 14.91 Denote the SCCSs passing through point Pij^k and projecting to 
C(e) as S{si),i = 1, . . . ,m. Assume that Q (Step 3 of Algorithm 14.91 ) is the set of isolation boxes 
of m points Qi, ■ ■ ■ ,Qs with Qi on S'(sj). Then in the plane x = 6j (or x = Oj), the surface becomes 
a plane curve f{bi,y, z) = and each surface patch passing through S{si) becomes a curve segment 
of the curve f{bi,y,z) = passing through Qi. We summarize this as the following lemma. 

Lemma 4.10 Use the above notations. If S is delineable over /(e) and I{c) respectively, then the 
number of TSPs passing through S{si) and projecting to R{c) is the number of curve branches 
passing through Qi and projecting to the region R{c). 

According to the above discussion, we have the following algorithm. 

Algorithm 4.11 NumTSP(/(x, y, z), Pij^k, e,c) S : f{x, y,z) = is the surface delineable over 
/(e) and I{c) respectively, Pij^k £ SV, e = {Pij, Pu,v) G cind c £ £C with e as an edge. The 
output is i^{Ei, c) where S{Ei) are all the SCCSes passing through Pij^k cmd projecting to C(e). 

1. Execute Algorithm NumSCCS(/(x, y, z),Pij^k, e). 

2. If e is not an x-vertical edge, execute the following steps 

(a) Let Q = {Qi, . . . , Qm} be the points obtained in Step 3 of Algorithm NumSCCS. 

(b) Let II21 = {g{bi,y), f{bi,y, z)} be the defining triangular system for Q. Execute Algo- 
rithm |33] with input Q to compute L#{Qi) and R^{Qi). 

(c) Let ci be the cell under e in the y direction and C2 the one above e. By Lemma 14.101 
#(5/, ci) = L#{Q^) and #{81,02) = R#{Qi). See Fig. [Elfor an illustration. 

3. If e is an x-vertical edge, execute the following steps 

(a) Let 7^ = {/?!,..., /?s} be the points obtained in Step 1 of Algorithm NumSCCS. 

(b) Let II22 = {dix^dij), f{x,dij, z)} be the defining triangular system for TZ. Execute 
Algorithm 13 . 3 1 with input TZ. 
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(c) Let ci be the cell on the left hand side of e and and C2 the cell on the right hand side of e. 
ByLemmaiJl #(r/,ci) = L#(i?/) and #(r,,C2) = R^Ri). 

If there exist no TSPs connect to a SCCS, then the SCCS is an isolated spatial curve segment. 
4.7 The General Case 

Until now, we assume that the surface S does not contain straight lines parallel to the z-axis. In this 
subsection, we will show how to treat surfaces that contain such lines. 

The aim is to get the points on the vertical lines where the topology of the surface changed, and the 
intersections between some SCCSes and the vertical lines, then the SCCSes originating from these 
points, and the surface patches originating from the line segments defined by these points. 

The following will show how to compute the special case when f{a,P,z) = for some point 
P = (a, (3). It is clear that g{x, y) = has a finite number of such points since f{x, y, z) has no 
factor containing x, y only. We can solve the problem in the following way. 

1. Take a coordinate system transformation such that the transformed line Li of the vertical line 
Lq can be projected as a line L2 on the new Xy-plane. 

2. Determine the topological information of L2: the intersections of L2 and the new projection 
curve, the number of curve segments originating from each intersection on its two sides. 

3. Determine the topological information of Li: lifting the intersections of L2 and the projection 
curve of the new surface to determine the corresponding points on Li. Find the points where 
the topology of surface changed on Li. 

4. We can made the same coordinate system transformation for the intersection of two surfaces 
G{x, y) = and f{x, y, z) = 0. Then we can decide the points on the vertical lines which are 
the intersections of SCCSes and the vertical line. 

5. Find the points where the topology of the original surface changed on Lq from Li by coordi- 
nate relationship. Determine the topological information of Lq. 

Remark: It is convenient to take a transformation such that L2 is a vertical line of the new projection 
curve if L2 can't overlap other line(s) of the projection curve. 

In this way, we can solve the special case in Algorithm 14.51 Since we have introduced the operations 
we need before, we just use an example to show the effectivity. 

In this special case, S8 contains the edge with the from {Pij^k^ Pi,j,k+i)- Similarly, SJ^ contain the 
face with the form {Pij,k, Pi,j,k+i,Pu,v,w)- 

We will continue the same example. Let us consider the following surface inside B = [—2,2] x 
[—2, 2] X [—2, 2] as an example. 

S : f{x, y, z) = x^y^ + x^z^ + y'^z^ - -xyz = 0. (16) 
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It is clear that /(O, 0, z) = 0. So (0, 0) is a point in special case. So Lq : {x = 0,y = 0, —2 < z < 
2}. The topology polyhedron of the surface is shown in Figure [141 

1. Take the system transformation 

{x = X - Z,y = Y - Z,z = Z.} (17) 

We get a new surface 

S':F{X,Y,Z) = X^Y^ -2X^YZ + 2X^Z^ -2XZY^ +4:XZ^Y - 4:XZ^ + 2Z^Y^ 
-4 Z^Y + 3 - I ZXY + I XZ^ + | YZ^ -\Z'^ 
= 0. 

Now Lo corresponds to the hne segment L\ : {X = Z ,Y = Z , —2 < Z < 2} on the new surface. 



Figure 14: Topology polyhedron 

Figure 13: Determining the in special case r r ^- i i- 

o t r ^ surface with vertical hne 

2. The projection curve of S' is shown in the right part of Figure \T3\ The red Une segment is the 
1/2 : {X — y = 0, — 2 < X < 2}. It corresponds to Li. The isolation boxes of the singularities of 
the projection curve of S' on L2 are below. 

[Pi,P2,P3,P,] := [[-|,-|]x[-|,-|],[-|i,-Sx[-i'-SaO,0]x[0,0],[|,|]x[|,|]]. 

3. The corresponding points Qi on Li of these singularities Pi, P2, Pi,, P^ are 

[<5l,Q2,Q3,Q4] := [[t X t X t], 

,7 7, , 41 655 , , , ,7 7,, 
' = [-4'-4l'l-i4'-1024l->°'°>->4'4ll' 
Assume the endpoints of Li are Qq, Q^. Computing the SCCSes originating from Qi{i = 0, . . . , 5) 
with Algorithm 14.91 we can find that: 



There are no SCCSes originating from QoiQb) except for QoQi {Q4Q5) on Li, we can find there are 
on surface patches originating from QoQi {Q4Q5) by Algorithm 14. 1 1 [ Similarly, Qi originates one 
sees from Li's two sides respectively, and the SCCSes originates two surface patches in the two 
cells besides Q1Q2, so Qi is a point we are interested; Q2 originates two SCCSes from Q2 besides 
Li respectively, and the four SCCSes all originate one surface patch on the cells besides Q1Q2 and 
Q2Q3, which means Q2 is not a point where the topology of the new surface changes on Li; Q3 
originate two line segments parallelhng to XY -plane as SCCSes on Li's two sides respectively, all 
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originating 2 surface patches on the cell bodies besides them, which means the SCCSes are singular 
curve of the surface, so Qs is a point we are interested; Q4 originate one SCCS on Li's two sides 
respectively, QsQ^ originates surface patches but there is no surface patches originating from 
so Q5 is a point we are interested. 

So we can conclude that Qi, Qs, are the points where the topology of the surface changed on Li. 

4. Take the same coordinate system transformation as ([TT] ) for g = xy{16x'^ + 16y^ — 49), We can 
get a surface: 

G{X,Y,Z) = {X - Z){Y - Z){16{X - Z)^ + IQ{Y - Zf - 49). 

We just need to decide some points on the line corresponding to the vertical line on the space curve 
defined by G{X, Y,Z) =0 and F{X, Y, Z) = 0. Use the method in II6I, we can find that there 
is only one point on the vertical line which is the intersection of SCCSes and the vertical line. It is 

[0,0] X [0,0] X [0,0]. 

5. Now we can get the points where we are interested on Lq, we can simply call these points as 
vertical points. By the coordinate relationship of Li and Lq, we can get the points we are interested 
on Lq which correspond to Qi. Since the topology of the surface does not change on Li at Q2, 
we need not to consider the corresponding point on Lq. Let Vq, Vi, V2, V3, V4 be the points on Lq 
corresponds to Qo,Qi,Q3,Q4, Q5 on Li. We have the points 

[^0,^l,^2,^3,V^4] = [[[0,0] X [0,0] Xt], 

t = [-2,-2],[-|,-|],[0,0],[|,|],[2,2]] 
and the edges {Vq, Vi), (14, ^2), {V2, V^), {V3, V^). 

Now we need to find out the SCCSes of the original surface which originate from these points and 
edges on vertical line. The basic idea is as below. 

At first, find a separate point Wi on each vertical edge, that is, between two adjacent points Vi-i, Vi, 
then construct a plane @Wi paralleling to XY-plane passing Wi, search a rectangle Ri containing Wi 
such that all the curve segments inside Ri originate from Wi, and when projected into XY -plane, all 
this kind of Ri correspond to a same rectangle R which only contains one critical point P. In order to 
determine the number of SCCSes originating from each vertical point, we need the following lemma. 

Lemma 4.12 The number of SCCSes originating from the point Vi equals the number of intersec- 
tions of line {x = a, y = /3} and the surface S between two planes <^f^d Qwi+x^ where {a, (3) is 
a point on C{e) inside R. 

Proof. Since there is only one vertical points between and @Wi+i^ the SCCSes between two 
planes originate from Vi. There is no part of the surface in Ri or has intersection with 1(e) 
when projected to R. Otherwise, there exists a critical point on R besides P. By Lemma l4~n the 
SCCS originating from Vi only intersects the line {x = a,y = (3} once. So the lemma is true. | 

So we have the method to decide the number of SCCSes for vertical line case. 

Then, we need to decide the number of surface patches originating from each vertical line segment 
in each plane cell. In fact, this is done! The boundaries of Ri have some intersections with the 
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surface, the number of the intersections in each cell body is the number of surface patches originat- 
ing from the corresponding vertical line segment. For our example, since Vq, V4 are endpoints and 
{Vq, Vi), (V3, V4) do not originate surface patch, we just need to find rectangles for (Vi, V2), (V2, V3). 
So we can conclude that {Vi , V2) originate two surface patches in cell bodies "2" and "4" respectively, 
and (V2, V3) originate two surface patches in cell bodies "1" and "3" respectively. When comput- 
ing the SCCSes originating from 14, V2, V3, we can find that Vi, V3 do not originate non- vertical 
SCCSes, V2 originates four line segments as SCCSes. 

In the end, we should form triangles for this case. The curve branches in Ri can intersect the plane 
triangles when projected to XY-p\a.ne. Use these points to subdivide the plane triangles, and then to 
form triangular patches. Note that when an endpoint of a plane triangle corresponding to a vertical 
line, some of the surface patches corresponding to the triangle should contain two or three TSPs. 

Figure [14] is a triangular polyhedron representation of the surface defined by Equation [16] which has 
a vertical line {x = 0,y = 0}. 

5 Ambient isotopic meshing of surface 

In this section, we will show how to compute an e-meshing of a surface S for a given e > 0. 

Let A^i be an e-meshing graph of the projection curve of S computed with Algorithm [STTOl Consider 
the two disjoint regions of B3 : 



Surface S has no singularities in the cylindrical region N3, so we can use a modified Pantinga-Vegter 
method [24] to compute its meshing. What we need to do is to compute the correct meshing inside 
S3. To present the algorithm, we need preparations given in Sections [5?T] and \52[ 

5.1 Extremal points of surfaces and spatial curves 

In order to give an ambient isotopic meshing for a surface, we need to consider 2:-extremal points of 
surfaces and spatial curves. A point is called z-extremal if the surface achieves a local extremum 
value at this point in the z-direction. We have 

Lemma 5.1 Let f{x, y, z) = ■ y, z) be a square free polynomial and fi irreducible polyno- 
mials. A necessary condition for the surface f{x, y,z) = to have a z-extremal point is 



S3 
N3 



B3 \ S3. 



(18) 
(19) 




(20) 



where only the nonzero resultants are included. 
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The following example shows that we need to consider the irreducible factors. Let / = {z — y){z — 
x){x'^ + + — 1). Then Res{f, fx, z) = Res(/, = 0. But the surface indeed has an 
z-extremal point at (0, 0, 1). 

Lemma 5.2 Let f{x,y,z) be a square free polynomial, D{x,y) defined in diOl ). Gi{x,y) defined 
in d20D . and r a fixed number Then D{x,y)G\{x,y) = is a necessary condition for the curves 
f{x,y,r) = 0, f{x,r,z) = 0, and f{r,y,z) = to have x-extremal, y-extremal, or z-extremal 
points. 

We also need to consider the z-extremal points of spatial curves defined by g{x, y) = f{x, y, z) = 0, 
where g and / are polynomials. For this purpose, we need to decompose the curve into irreducible 
ones. The leading coefficient of g (/) as an univariate polynomial in y {z) is called the initial of g 
(/). Two polynomials of the form g{x, y), f{x, y, z) is called an irreducible chain if the following 
conditions are satisfied | |2T1 (pages, 297-381) 

• g{x, y) is an irreducible polynomial. 

• f{x,y,z) is an irreducible polynomial of z module g = 0, deg(/, y) < deg{g,y), and the 
initial of / is a polynomial in x. 

For instance, g = y"^ — x, f = z"^ — x is not irreducible, since / = {z — y){z + y)+g = {z — y){z + y) 
mod {g). 

For an irreducible chain g{x, y), f{x, y, z), we define its saturation ideal to be 

Sat(g(x, y), f{x, y, z)) = {P \ Il^P G (/, g)} 

where Ii and I2 aie the initials of g and / respectively. It is known that the saturation ideal of an 
irreducible chain is a prime ideal, and thus defines an irreducible spatial curve 1211 (pages, 297-381). 

Any spatial curve f{x, y, z) = g{x, y) = can be decomposed into the union of irreducible curves 
algorithmically: 

V{g{x, y), fix, y, z)) = U,F(Sat(y,(x, y), f\{x, y, z))) (21) 
where gi{x, y), fi{x, y, z) are irreducible chains. We can prove the following result: 
Theorem 5.3 Let g{x, y), f{x, y, z) be an irreducible chain and 

I{x) = product of the initials of f, g. (22) 
T{x) Res{Res{h,f,z),g,y) where h{x,y,z) ^ f^gy- fyg^. 

Let E be the set of z-extremal points of the curve C : f = g = 0. Then 

ProUE) C V{T{x)) U V{I{x)). (23) 

Furthermore, ifT{x) = 0, then the curve y(Sat({/, g})) is contained in several planes perpendicu- 
lar to the z-axis. 

Proof. For any point P = {a, [5, 7) on C, the necessary condition of P being a z-extremal point of 
C is the tangent line of C at P is perpendicular to z-axes. If P is neither the singular point of / = 
nor (7 = 0, the tangent planes of / = and (7 = at P are both well defined. The tangent line of C 
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at P is the intersection of the tangent planes of / = and g = at P. The tangent direction n of C 
at Pis 

n = <fx,fy,fz>\pX<gx,gy,0>\p 
= < -9yfz, fzOx, fxQy - fy9x > \p- 

Since the tangent planes of / = and (7 = at P are: 

( (x- a)Ma, 7) + (y - (a, A 7) + - 7)/.(a, l) = 0, 
\ {x-a)ga;{a,(3) + {y- f3)gy{a,f3)=0 

respectively, n is perpendicular to the z-axes, that is: 

n-<0,0, 1> = f,^{a,P,'y)gy{a,(3) - fy{a,P,-f)g^{a,(3) 
= h{a,(3,-f) =0. 

Therefore, E C V{h{x, y, z)). (|23]l is true. 

If T{x) = 0, we prove that V{{f, g} / 1) is contained in several planes perpendicular to the z-axis. 

First we claim that n is well defined on C except finite number of points, that is only finite number 
of points on C are the singular points of / = or 51 = 0. If it is not true, at least one of the following 
conditions occurs: 

CI. V{f, g, fx, fy, fz) has 1-dimensional component. 

C2. V{f, g, gx,gy) has 1-dimensional component. 

If CI. occurs, it means that fz G Sat(/, g). It is impossible. Condition C2. could not take place for 
the same reason. Note that h{xo,yo, zq) = for any point {xo,yo, zq) on C. The tangent direction of 
C at almost all points is the form {A, B,0). 

Then we prove this component of C lies in some planes z = zq. This component of C can be 
parametrization in some segments. Assume the parametric equation is r(t). We have 

r(t) =r(to)+ / r'{t) = {x{t),y{t), zq), 
Jo 

where r(to) = {xQ,yQ, zq). It implies that this segment of C lies in the plane z = zq. Therefore, 
the irreducible component of C which contains this segment lies in the plane z = zq. We prove this 
theorem. I 

The following example shows that we need to decompose the curve into irreducible ones. Let / = 
z{x'^ + z"^ — I), g = y. Then Res{fxgy — fygx, f, z) = 0. But the curve indeed has a z-extremal 
point at (0,0, 1). 

5.2 Compute segregating box for an SCCS 

In Section 1431 we showed how to compute the segregating box for a singular point. In this section, 
we introduce the concept of segregating boxes for singular curve segments. 
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In Algorithm 13.101 a curve segment C(e) of the projection curve C is represented by a segment e 
contained in a box Be, as shown in Fig. |5] When lifting C to the space, we obtain a set of SCCSes 
Si, i = 1, ... ,d of S represented by edges Ei G S£ (see Section 1421 ). A box Bg. = Bg x [cj, fi] is 
called a segregating box for Si if B^^ n B^^ =0 for i ^ j and S does not intersect with the top and 
bottom faces of B5.. In Fig. \T5\ we give a segregating box for the surface patches AiBiCiDi and 
A2B2C2D2 intersecting at curve segment P11P21 which is lifted from curve segment PiP2- 

Assume that all Si are monotonous in the direction of z, the following algorithm shows how to 
compute segregating boxes for the SCCSes: Si,i = 1, . . . ,d. 





Figure 15: Segregating boxes 



Figure 16: Merge meshes 



Figure 17: Divide N2 into 
boxes 



Algorithm 5.4 SegBoxC(/(x, y, z), g(a;, y), B3, Bg, e). Let S : f{x,y,z) = be the surface, B3 
the bounding box, Bg a nice box (see Fig. E]) containing a curve segment C{e) of the projection curve 
C:g{x,y)=OofS,e>0. 

The output is a pair (P, S). P is a set of interior-disjoint boxes contained in Bg, the union of which 
contains C(e). For each P j G P, there exist 3D boxes Sij G S which are the segregating boxes for 
the SCCSes lifted from C(e) H P^. 

1. We consider case (a) in Fig. IH Other cases can be treated similarly. C(e) divides Bg into two 
cells ci and 02- 

2. Let Pi = {Bg}, P = 0, S = 0. Repeat the following steps until Pi = 0. 

(a) Let B = [a, b] x [c, d] G Pi and remove B from Pi. 

(b) Execute RootIsol({(7(a, y), f{a, y,z)}, [c, d] x [Zi , Z2], e) to compute the points Pi^i,i = 
1, . . . , A^i lifted from Pi. See Fig. [T5]for an illustration. Let the isolation box for Pi j be 

Si,i X [ei,j, /i,j]. 

(c) Similarly, let P2,i; 62,4, /2,i, ^ = 1, • • • , be the points Ufted from P2. By Lemma 1411 

Ni = N2. 

(d) LetBi = BgX [min{ei,j, 62,*}, max{/i,i, /2,i}], i = 1, . . . , A^i. 

(e) If |Bi| < e for all i, add Bg to P and add Bj to S. 

(f) Otherwise, subdivide Bg into four equal boxes and add the boxes intersecting with C into 
Pi- 



3. Repeat the following steps until all boxes in S are segregating. 
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(a) Let B = B X [e, /] € S. 

(b) If 0/(6 X [e, e]) and □/(B x [/, /]), then B is segregating, we do nothing for 
B. 

(c) Otherwise, remove B from P and B from S. Subdivide B into four equal boxes Ci, C2, 
C3, C4, add each C, intersecting with C into P, and add Cj x [e, /] into S. 

4. Return P and S. 

Proof of correctness. We need only prove the termination of the algorithm. According to the as- 
sumption, all Si are monotonous in the z direction. So Step 2 terminates in a finite number of steps. 

At the beginning of Step 3, for any C(e) C B = [a,b] x [c,d] £ P where e = (^1,^2), -Pi = 
(a, a), i-*2 = {b, (3), let C{sj) be a curve segment lifted from e and Bj = [o, b] x [c, d] x [cj, fi] be 
the corresponding box where Sj = {Pi,i,P2,i), Pi,i = {a, a, 7^), P2,i = P, Ti), we have |Bj| < e 
and 

/(a, a, ei)/(a, a, ei)/(^ /?, fi) / 0. 

Furthermore, does not intersect with the top nor bottom faces of Bj since Si is monotonous in the 
direction of z. That is □/(C(e) x [ej, ei])f{C{e) x [/j, /j]). So there exists a positive number 5 
such that 

^□/(d(C(e),<5) X [e,,e,])/(d(C(e),<5) x [/^„/,]) 

where (i(C(e), 5) is a zonal region in containing the points Q such that the distance between Q 
and C(e) is less than 5. We can get the set of sub-boxes of B in a finite steps such that all boxes in it 
are contained in the region d{C{e), 5). Then the algorithm clearly terminates. | 

5.3 Compute e-meshing of surface 

Similar to the case of curves, we need to modify the Pantinga-Vegter method. A box B is called a 
nice box if each face of B is a nice 2D box. For an illustration, see the 2D case in Fig. ID To make 
the process precise, we introduce the following definition. 

A meshing polyliedron of a surface 5 is a four-tuple J\A = {V, £, T ^ IS] where ("P, <S, T) is a 
polyhedron whose vertices are with rational numbers as coordinates and whose faces are the meshes 
for 5; ;B is a set of nice boxes and segregating boxes of singular points of S s.t. for each F £ T, 
there exists a Bi? G ^ with the property that the surface patch 5 H Bi? is connected. 

A meshing polyhedron M. is called an e-meshing polyhedron if each box B in B satisfies |B| < e. It 
is easy to show that an e-meshing polyhedron for a surface S provides an e-meshing for S according 
to the definition given in Section 2. 

Algorithm 5.5 MPV3(/(x, y, z), N3, e). S : f{x, y,z) =0 is the surface. N3 is a box contains no 
zero ofD{x, y)Gi (x, y) = where D{x,y) is defined in diOD and Gi {x, y) is defined in d20D . Output 
an e-meshing polyhedron for ST<s,y 

1 . Subdivide N3 into boxes B j at the comer lines (Fig. [17] shows how to subdivide the region 
N2 in Fig. Oa), where the dotted lines are newly added.) and execute the Pantinga-Vegter 
algorithm with initial values {Bj}. Let S be the output. 
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2. For each cube B G S, repeat subdividing B until all of the following statements are false. 

(a) There exists an edge {A, B) of B s.t. G □/((A, B)) and f{A)f{B) > 0. 

(b) There exists a face ABCD of B s.t. f{A)f{B) < A f{B)f{C) < A f{C)f{D) < 
OA f{D) f {A) <0. 

(c) |B| > e. 

Termination of the algorithm is guaranteed by Lemma [5^ 
Now we can compute the e-meshing for . 

Algorithm 5.6 ATopSur(/(x, y, z), B3, e). The input is the same as Alsorithm \4.5\ The output is 
an e-meshing polyhedron /or ^Bs. 

51 Compute the critical points of the projection curve and their segregating boxes. 

1. Let 

G{x,y) = sqrfree(G(x,y)Gi(2;,y)), (24) 
where G is defined in (fTTl) and Gi is defined in (l20l ). 

2. Execute the first four steps of Algorithm 13.51 with input (G(x, y), B2, e) to compute a 
set of points Vi and the segregating box for each point in Vi. We need to modify the 
algorithm as follows. In Step 3 of Algorithm l3.5[ we use the new projection polynomial: 

H{x) := H{x) n Ii{x) n Ti{x), (25) 

i i 

where H{x) is defined in (|4|, Ii{x) and Tj(x) are defined in (l22l ) with decomposition 
((2TI) . Only the nonzero Tj are considered. 

52 Compute SVq and the set SBq of segregating boxes for points in SVq- For any Pij e Vi, 

use Algorithm 14.71 with input (/, B3, j, e) to compute the points lifted from Bij and their 
segregating boxes. Let Bi be the set of all updated segregating boxes Bjj- of Pij. Let Aii = 

{ri,Bi}. 

53 Compute an e-meshing graph for the non-singular part of in N2 defined in dill. Let 

Mo = MPV2(G(x, y), N2, e), where N2 is defined in ©. 

54 Compute segregating boxes for SCCSs: 

1. Assume Mq = {Vo,£o,Bo}. Let SBi =V2 = £2 = B2 = 9. 

2. For each B € Bq, execute the following steps: 

(a) Compute {P,S}=SegBoxC(/(x,y,z),G(x,y), B3,B,e)|!] 

(b) 5^1 = 5^1 U S and update V2,£2, B2 according to P which subdivides B. 

3. LetM2 = {V2,£2,B2}. 

55 Compute the extended meshing graph £Gs of C with Algorithm ?? with input Mi and M2- 

'step SI ensures all z-extremal points of the curve C : f = 0, G — are in S3. Hence the SCCS in B is monotonous 
in the direction of z. 
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S6 Meshing the singular part of 5 in S3. 



1. Let {S'Pi,S£i,SJ'i}=TopSur{f{x, y, z), B3). Modify Algorithm TopSur as follows: 
use G{x, y) defined in ((24l) in Step 1, use £Qs in the Step 2, and use SVq in Step 3. We 
actually only run Steps 4 and 5 of Algorithm TopSur. 

2. Let TWi = {SVi^SSi^SJ'i^SBq U SBi] where 5^0 and SBi are from Steps S2 and 
S4 respectively. 1 is an e-meshing polyhedron for . 

57 Meshing the non-singular part of S in N3. Let M2 = {MV2, M£2, MT2, MB2} = 

MPV3(/(x, y, z), N3, e), where N3 is defined in 

58 Merge A^i and to obtain an e-meshing polyhedron for S. Output Merge (A^ 1 , M2) (with 

Algorithm EH). 

Theorem 5.7 Algorithm 15. 61 computes an e-AIMESH for Sbs- 

Proof. The prove is similar to the proof of Theorem l3.11l | 

In principle, there exist no difficulties to implement the algorithm. But, it will take a lots of time, since 
we need to incorporate algorithms from symbolic computation, interval arithmetics, and marching 
cube into one program. This will be our further work. 

In the final step of Algorithm 15. 6[ we need to merge two meshing polyhedrons, which will be done 
by the following algorithm. 

Algorithm 5.8 Merge(A4i, A^2)- Mi = {MVi,M£i,MJ'i, MBi} and M2 = {MV2,M£2, 
A4J-2,MB2} are the e-meshing polyhedrons ofS in S3 and N3 respectively. The algorithm merges 
^Al and M2 ond outputs an e-meshing polyhedron M = {MV, MS, MJ-} for the surface. 

51 Let MBt = MBi. 

52 While MBt + 0, repeat 

1. Remove B = [a, 6] x [c, (i] x [e,/] from MBt. Insert box Bj = [6, 6i] x [cj,dj] x 

G MB2 which is connected with B according to S and adjacent to the face 
F = [fe, 6] X [c, d] X [e, /] into Bq and insert corresponding ^(Bj) into Va ■ Pick out 
boxes Bj satisfying di = miriB eBai'^i}- Rename them to be Bi, . . . ,Bm. Sort the 
residual boxes in B^ as {Bm+i, • • • , B^} such that Cm+i < Cm+2 < ■ ■ ■ < Cr and for 
each Bfc, k > m, B^ is connected with some Bj, j < m according to 5(Note that the 
result is not unique, and any Bk,0 < k < m only overlaps with B on the vertical edge 
[b,b] X [di,di] X [ciji]). 

2. For i from 1 to r do 

(a) Remove the points P from V{Bi) and insert points Q £ SV n (B n Bj) into T^(Bj) 
if P G L U R where L = [6, h\ x [c, c] x [e, /], R = [6, b] x [d, d\ x [e, /]. Remove 
edge (P, Pi) from M£2 which are the edges with P as an ending point and insert 
{Q, Pi) into MJ^2- Remove triangular faces (P, Pi,Pj) from MJ^2 which are the 
faces with {P,Pi) as an edge and insert {Q,Pi,Pj) into MJ^2- 
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(b) If there exist P G Ri where Ri = [b, b] x [di,di] x [e^, Cj] and Ri n L n R = 0, add 
Pinto MVi. goto (e). 

(c) If there existP G Di C F where Di = [b, b] x [q, dj] x [e^, Cj]. Add P in 5^. goto 
(e). 

(d) If there existP G Ui C F where Ui = [b, b] x [ci,di] x [fi, fi]. Add P in cSP. goto 
(e). 

(e) Assume the other point contained in the face [6, b] x [cj, di] x [cj, /«] of Bj is Q and 
(Q, S, T) G TWJTi is the triangular face with Q as a vertex where 5 G R. Remove 
{Q,S) from M.£i and insert {Q,P), {P, S) into MSi. Remove triangular faces 
(Q, 5, T) from MJ='i and insert (Q, P, 5), (P, 5, T) into J"i(The four edges of 
this face of Bj contains two points. We can always assume that we have dealt with 
the other one point, since Bj, z > m is connected with some B^ we have dealt with). 

(f) Update M.B2 according to the new F(Bj). 

3. Determine the connection information of the other three faces of B in the similar way. 

S3 Out M = {MV, MS, MJ^} where MV = MVi U MV2, MS = MSi U MS2, and MJ^ = 
MJ^i U MF2- 

A box Bi is said to be adjacent to a box B2 w.r.t. the surface S if Bi and B2 are interiorly disjoint 
and S intersects Bi n B2. We need only consider how to merge the meshes in two adjacent boxes. 

We use the example in Fig. [16] to explain the algorithm. The large box B is in S3 and contains 
singularities. We consider the right face F of B. Let Bj, i = 1, . . . , 5 be the boxes in N3 adjacent to 
B at face F. By Step 1 of Algorithm l5.5[ Bj must be completely between lines AB and CD. We will 
adjust the meshes in B and leave the meshes in Bj unchanged. Since all the meshes are triangular, 
let OPQ be the mesh of S in B, and NiPi^iPi the mesh of S in Bj. We will replace the mesh OPQ 
with the meshes Afj = OPi-iPi,i = 0, . . . , 4. If Pj is above BC, Pi is taken to be the intersection 
of BC and the line passing through Pj and parallel to AB. Other cases can be treated similarly. 

6 Conclusion 

This paper proposes complete methods to compute isotopic and ambient isotopic meshings for im- 
plicit algebraic curves and surfaces. We use symbolic computation to achieve completeness and 
whenever possible use interval arithmetics to achieve practical effectiveness. Note that an isotopic 
meshing without precision and an e-meshing are quite different and can be used for different pur- 
poses. 
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